One dimensional heat conductivity exponent from random collision model 
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We study numerically the thermal conductivity coefficient n as a function of system length L for 
several different quasi one dimensional models: classical gases of hard spheres with both longitudinal 
and transverse degrees of freedom. We introduce a model that is ergodic and highly chaotic but 
also conserves energy and momentum, and is very useful because it shows scaling even at small 
system sizes. We find that k ~ L a over more than two decades, with a very close to the analytical 
prediction of 1/3. 
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Since the surprising result obtained over thirty 
years ago that the heat current flowing across a one- 
dimensional chain of harmonic oscillators with a small 
temperature difference between the two ends is indepen- 
dent of the length L of the chain Q], the conductivity 
of one dimensional systems has been studied analyti- 
cally dHli and numerically 0, H H HQ at great 
length. The standard approach to conductivity would 
predict that if the temperature gradient VT in a mate- 
rial is small, the heat current flowing through should be 
of the form j = -kVT, where k is a property of the ma- 
terial. On the other hand, the result for the harmonic 
oscillator chain, that a small temperature difference AT 
results in a current j oc AT that is independent of L, is 
equivalent to a conductivity k oc L. In the various models 
studied thereafter, singular conductivities k ~ L a have 
been found, with a variety of possible values for a. On 
the other hand, for some one-dimensional models a con- 
ventional L-independent k has also been obtained. 

It is now believed that the singular conductivity of 
these models has two possible causes. Firstly, if the 
model is integrable, as in the case of the harmonic os- 
cillator chain, the system does not equilibrate thermally. 
The behavior of the conductivity then depends on the 
details of the system. In fact, it has been shown that by 
changing the coupling of the oscillator chain to the heat 
reservoirs at the ends, normally a benign procedure, one 
can tune the exponent a over a range [3| ■ Secondly, even 
if the model is not integrable, if the internal interactions 
in the system conserve momentum, the conductivity is 
singular, due to advection of heat in long wavelength 
modes 0, H 0] . If a model is not integrable and does 
not conserve momentum, k should have a well defined 
limit as L diverges ^1 • Analytical studies 0, 0] predict 
that for non-integrable momentum conserving systems, 
a should have a universal value of 1/3. 

Numerical results for momentum conserving systems 
have yielded values of a ranging from 0.25 to 0.5. Recent 
studies of one dimensional chains of hard point particles 
with alternating masses have shown unexpectedly large 
corrections to scaling even for systems of ~ 10 4 particles, 
with a estimated to be 0.25 Q and 0.33 p. Similar re- 
sults have been obtained for chains of Fermi Pasta Ulam 
chains, with a estimated to be 0.37 For the system of 



hard point particles, the slow convergence to the asymp- 
totic behavior has been justified [6| by noticing that the 
system is always 'close' to an integrable model. Thus 
if one considers a chain of particles with equal masses, 
energy is carried ballistically, the system does not ther- 
malize, and a — I. On the other hand, if the ratio of the 
masses of successive particles in the alternating chain is 
chosen to be very different from unity, the light particles 
are almost inconsequential, and the problem again re- 
duces to one of equal mass particles. In Ref. 6], a mass 
ratio of 2.62 was found to yield the longest power-law 
scaling range, from which a was estimated. 

In this paper, we study how to eliminate residual ef- 
fects of integrability by considering non-integrable and 
highly chaotic models. We start with Sinai and Chernov's 
pencase model that was initially introduced by them 
to study one dimensional hydrodynamics. In this model, 
hard sphere particles are confined to a long narrow tube. 
We consider both periodic and hard wall boundary con- 
ditions in the transverse direction, and apply heat baths 
to the two ends. The extent in the transverse direction 
is taken to be slightly less than twice the diameter of the 
particles. This ensures that the particles cannot get past 
each other, but allows a large range of incidence angles 
at the collisions. Thus the transport of energy along the 
tube remains quasi one dimensional, with the transverse 
degree of freedom serving as an additional randomizing 
effect. This model (without heat reservoirs) has been 
proved to be ergodic in four dimensions and hyperbolic 
in three . 

The extra degree of freedom in our model limits the 
system sizes we simulate. For the range of system sizes 
that we consider, there is still insufficient universality 
to obtain the asymptotic value for a with confidence. 
Therefore, we have devised a novel Id model that we call 
the "random collision" (RC) model, discussed in detail 
below, that is extremely chaotic, yet still satisfies energy 
and momentum conservation. 

Probably due to the further randomization introduced 
by our new dynamics, the conductivity fits very well to 
the form k ~ L a , with a close to 1/3. If the masses of all 
the particles are taken to be equal, the estimated value of 
a is 0.29 ±0.01, while when the particle masses alternate 
with a mass ratio of 2.62, one obtains a to be 0.335±0.01. 
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The small discrepancy from the theoretical prediction of 
a = 1/3, although larger than the error bars, is within 
the range one might expect from corrections to scaling 
from the leading irrelevant operators. 

We also show results for the autocorrelation function 
of the energy current, which is related to the conductiv- 
ity through the Kubo formula ^4|. The autocorrelation 
function has a zero frequency limit that scales close to y 
L 1 / 3 , in accordance with previous predictions and our 
numerical results for the conductivity. However, the full 
correlation function C(oj, L) shows quite complicated be- 
havior, not easily described by any simple scaling form, 
possibly indicating the existence of multiple time scales. 

One might be concerned whether the analytical deriva- 
tion of a = 1/3, which uses the hydrodynamic descrip- 
tion of a one-dimensional normal fluid [lOj is valid for 
the models we have considered here. Since the trans- 
verse direction is small, it should not affect long wave- 
length singularities in the dynamics. Of greater concern 
is, with periodic boundary conditions, the existence of 
the transverse momentum as another conserved quantity. 
Although this makes the hydrodynamic equations more 
complicated, and increases their number from three to 
four, we recall that the analytical calculation relies only 
on the fact that (without an applied temperature gradi- 
ent) the system reaches thermal equilibrium, that it sat- 
isfies Galilean invariance, and that a finite sound velocity 
sets a cutoff to the dynamics for any finite system size. 
None of these conditions is violated by the introduction 
of the transverse momentum. 

Figure ^ shows a log log plot of the conductivity as 
a function of the length the system L, for the pencase 
model. The results for both periodic and hard wall 
boundary conditions are shown. All the particles were 
taken to have the same mass. The diameter of the parti- 
cles was 0.6, while the cross section of the tube and the 
average longitudinal separation between centers of neigh- 
boring particles were both taken to be 1. The temper- 
atures of the heat reservoirs at the two ends were taken 
to be 1.0 and 1.2 respectively; it was verified numeri- 
cally that this temperature difference is sufficiently small 
for the system to be in the linear response regime. The 
heat reservoirs at the ends were implemented as follows: 
whenever an extremal particle collided with the reservoir 
adjoining it, its velocity was randomized, drawn from the 
distribution P(v x ,v y ) oc v x exp[— m(v x + Vy)/(2kBT)] 1 
where x and y are along the longitudinal and transverse 
direction respectively and T is the temperature of the 
reservoir. (This is the velocity distribution for particles 
leaking out of a heat reservoir.) At each such collision, 
the energy exchanged by the system and the reservoir is 
kept track of, and used to calculate the time average of 
the energy current at both ends. 

In the same figure, the conductivity as a function of 
system length is also shown for a different choice of model 
parameters: alternating particle masses, with a mass ra- 
tio of 2.62, a tube cross section of 1.14 and a longitudinal 
interparticle separation of 0.9. (Only hard wall trans- 
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FIG. 1: Log-log plot of the conductivity as a function of the 
number of particles for Sinai's pencase model I I J . The upper 
two plots are for periodic and hard wall boundary conditions 
in the transverse direction, with the mass of all particles be- 
ing 1. The slopes of these plots, which should be equal to the 
conductivity exponent a, are 0.25±0.02 and 0.26±0.01 respec- 
tively. The lowest plot is for hard wall boundary conditions in 
the transverse direction, with the mass of particles alternat- 
ing between 2.62 and 1. The slope of the plot is 0.34 ± 0.02. 
Even though each individual plot fits well to a straight line 
(apart for deviations at the low L end), the differences be- 
tween the slopes preclude a good estimate of the asymptotic 
large L value of a. 

verse boundary conditions are shown for this case.) 

The number of particles in the system ranged from 8 
to 1024. This is substantially less than the largest sys- 
tem sizes used when simulating the one dimensional chain 
of hard point particles. However, introducing the trans- 
verse dimension should make the system no longer near 
integrability, and therefore allow the large L limit to be 
reached quickly. Unfortunately, as seen in Figure Q this 
is not the case. There is some curvature in all the plots; 
more importantly, there is substantial disagreement be- 
tween the slopes obtained when the particles have equal 
or alternating masses. Note that the plots all curve down- 
wards, from which one might be tempted to conclude 
that the asymptotic slope (i.e. the value of a) would be 
smaller than obtained from the curves. However, prior 
experience with the one dimensional system [(| indicates 
that it is possible for the curves to turn around at much 
larger system sizes. Thus one cannot obtain even an up- 
per bound to a from the figure, and must conclude that 
the corrections to scaling are large for this model |l5|. 

In order to randomize the dynamics further, enabling 
faster convergence to the asymptotic scaling form for 
large L, we modify the model above. Firstly, the di- 
ameter of the particles is taken to be negligibly small, 
while keeping the cross-section of the tube as less than 
twice the diameter. Secondly, the particles are no longer 
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FIG. 2: Log-log plot of the conductivity as a function of 
the number of particles for the RC model introduced in this 
paper. The upper plot is for all the particles with mass 1, 
while the lower plot is for particles whose masses alternate 
between 1 and 2.62. The slopes for the two plots are 0.29±0.01 
and 0.335 ± 0.01. Both the plots are for periodic boundary 
conditions in the transverse direction; the results for hard 
wall boundary conditions were very similar. 



disk shaped, but highly irregular. As a result, when two 
particles collide with each other, in the center of mass 
frame they can recoil in any direction, unrelated to the 
direction of impact. For any collision, we take the recoil 
angle in the center of mass frame to be a uniform random 
variable ^(| . The result of both these modifications to- 
gether is that the transverse coordinate y of the particles 
becomes redundant, and they effectively move along the 
x axis. However, the transverse velocities v y are retained. 
In this RC model, each particle has both v x and v y , with 
the latter behaving as an auxiliary variable that is only 
important in collisions. In any interparticlc collision, the 
total momentum in the x and y directions and the to- 
tal energy are conserved. Collisions with the reservoirs 
at the two ends are still implemented as before. In the 
transverse direction, periodic boundary conditions cor- 
respond to v y for a particle remaining constant between 
collisions, whereas hard wall boundary conditions allow 
v y to be reversed. In the latter case, since in the zero 
cross-section limit any particle undergoes a huge number 
of collisions with the sidewalls between two collisions with 
its neighbors, one should change the sign of v y randomly 
between interparticle collisions. 

Figure [3 shows a log log plot of the conductivity as a 
function of L for the RC model. The particles are now 
points, the average interparticle separation is unity, and 
the temperatures of the reservoirs are 1 and 1.2 as be- 
fore. In the linear response regime, the dependence of 
the energy current on the average interparticle separa- 
tion and the reservoir temperatures can be found triv- 
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FIG. 3: Scaling plot for the autocorrelation function of the 
energy current, C(io), as a function of frequency u) and system 
size L. The x-axis is ujL and the y axis is L~ 1//3 C(a;, L), 
in accordance with analytical predictions. All particles have 
equal mass, and periodic boundary conditions are imposed 
in the longitudinal direction. The low frequency limit of the 
scaled functions approaches an L-independent constant. 



ially, so the only parameters that can be varied (apart 
from the number of particles) are the masses of the par- 
ticles. Only the results for periodic boundary conditions 
in the transverse direction are shown, as the results for 
the hard-wall case are very similar. As before, plots are 
shown for the case when all particles have the same mass 
and the case when the particle masses alternate with a 
mass ratio of 2.62. Both parameter choices yield plots 
that rapidly converge to a power law form, but with 
slightly different exponents. The estimated exponents 
are a = 0.29 ± 0.01 and a = 0.335 ± 0.01 respectively 
for the two cases. Both of these are very close to the 
analytical prediction of a — 1/3. Although the difference 
between the two estimates for a is larger than the error 
bars, we expect that this is due to corrections to scal- 
ing from irrelevant operators in a renormalization group 
analysis; an O(l) bare strength for irrelevant operators 
can produce effective values of a that differ from 1/3 
by the desired amount. Thus the numerical results of 
Figure |3 are a strong indication of the validity of the 
prediction of a = 1/3. 

We have also measured the autocorrelation function 
for the energy current, C(t) = (J(t + t')J(t')), where the 
average is over t', since this is related to the conductivity 
by the Kubo formula 0]. This was done in the center 
of mass frame, as necessary for the Kubo formula |17| . 
Figure [3] shows the Fourier transform of the correlation 
function, C{u>), with periodic boundary conditions in the 
x direction (without heat reservoirs), for the case when 
all the particles have identical mass. Based on the ana- 
lytical prediction [T^j, we show L~ 1 l z C(u>L) for a range 
of system sizes. The u> — > limit of all the curves ap- 
pears to be the same, in accordance with the numerical 
results we have obtained for the conductivity, although 
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FIG. 4: Unsealed plot of the autocorrelation function of the 
energy current, C(w), as a function of frequency uj. The upper 
curve corresponds to a system with 512 particles and the lower 
curve to a system of 1024 particles. The two curves have been 
shifted with respect to each other for clarity. All particles have 
equal mass, and periodic boundary conditions are imposed in 
the longitudinal direction. The straight lines correspond to a 
~ a; -1 / 4 functional form. 

the error bars are clearly larger due to the method of 
measurement. The collapse of the curves to a single plot 
at low frequencies is less clear; in fact, as shown in Fig- 
ure 01 a functional form of C(lu,L) ~ w -1 / 4 seems to 
work slightly better. The correlation function flattens 
out from this form at low frequencies, at a cutoff ujl that 
appears (for L ranging from 64 to 1024, not all shown in 
the figure) to decrease slightly faster than ~ 1/L with L. 
The existence of a well defined limit L _1 / 3 C(o> — > 0,L) 



would require a 1/L 4 / 3 decay for this cutoff scale. In ad- 
dition to this cutoff, there is a clear resonance peak at 
lj ~ 1/L, indicating that there are important dynamics 
on the t ~ L timescale. 

It is not clear whether the correlation function con- 
verges to the simple scaling form for much larger system 
sizes, or whether there are indeed multiple time scales in 
the dynamics: t\{L) ~ L and t2(L) ~ L 13 with j3 rs 4/3. 
However, earlier work 0] on the one dimensional KPZ 
equation [Tflj and work in two dimensions pfij , has found 
it difficult to obtain universal critical exponents from nu- 
merical simulations, even with larg e system sizes. The 
hydrodynamic description used [lfj to obtain a is simi- 
lar to the KPZ (Burgers) equation, but with three equa- 
tions instead of one. In fact, a was correctly estimated 
earlier _(| from the KPZ equation. Thus the slow conver- 
gence that we see for a for the pencase model — and, to a 
lesser extent, for the RC model — may be a similar phe- 
nomenon to that seen for the KPZ equation. Of course, 
even the hydrodynamic description for one-dimensional 
conduction does involve two time scales: for propa- 
gation and decay of the sound modes. 

In this paper, we have introduced a random collision 
model for studying dynamics of momentum conserving 
one dimensional systems, in order to obtain the scaling 
form of the thermal conductivity function of system 

size L. This has allowed us to probe the scaling regime 
better than has been possible with previous models that 
show a slower approach to the large L limit and less ro- 
bust behavior. Over a wide range of length scales, we find 
good agreement with the earlier analytical prediction of 

We thank Abhishek Dhar and M.A. Moore for very 
useful discussions. 
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